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ABSTRACT 

We continue our earlier work on a class on nonlinearly stable Runge-Kutta local projec- 
tion discontinuous Galerkin (RKDG) finite element methods for conservation laws. Two- 
dimensional Euler equations for gas dynamics are solved using P 1 elements. We discuss 
the generalization of the local projection, which for scalar nonlinear conservation laws was 
designed to satisfy a local maximum principle, to systems of conservation laws such as the 
Euler equations of gas dynamics using local characteristic decompositions. Numerical exam- 
ples include the standard regular shock reflection problem, the forward facing step problem 
and the double Mach reflection problem. These preliminary numerical examples are chosen 
to show the capacity of our approach to obtain nonlinearly stable results comparable with 
the modern nonoscillatory finite difference methods. Generalizations to P* elements with 
k > 1 and the use of adaptive triangulations to minimize local errors constitute ongoing 
research. 
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1. Introduction. In this paper we continue our earlier work [Cl], [C2], [C3], [C4], 
of constructing and analyzing a class of discontinuous Galerkin finite element method for 
solving conservation laws 

d 

(1.1) u t + 52(f;(u)) Xi . = °> 

i=l 

equipped with suitable initial or initial-boundary conditions. We concentrate on two- 
dimensional Euler equations of gas dynamics, i.e., in (1.1) with d = 2, u = (p, pq z , pq y , E) 1 , 
fj (u) = q z u + (0,p,0,q x py, f 2 (u) = ^ y u-|-(0,0,p,g y p) t , where q x ,q y axe the velocity com- 
ponents in the x and y directions, p is the density, E is the total energy, 
p = (7 — l)(E — \p{q 2 x + <Z y )) is the pressure, and 7 = 1.4 for the air. 

One distinctive feature of our approach is a local projection limiting which borrows the 
successful non-oscillatory finite difference methodology, guarantees total variation bound- 
edness (TVB) for one- dimensional nonlinear scalar equations and linear systems, and yields 
a local maximum principle for multi-dimensional scalar equations. Another feature of our 
approach is the use of high-order total variation diminishing (TVD) Runge-Kutta type 
time discretizations [Si] which renders the scheme explicit (and hence fully parallelizable) 
and computationally efficient. The general framework was. given in [C2] for the nonlinear 
one-dimensional scalar case. The TVB property and convergence were proven for general 
P* elements (which give rise to uniformly ( k + l)-th order accurate schemes). Numeri- 
cal results for k = 1 and k = 2 (second and third order) were shown which gave sharp, 
monotone shock transitions and uniform high-order accuracy in the smooth part of the 
solutions. In [C3], we applied the method to the one-dimensional Euler equations of gas 
dynamics by designing the local projection in the local characteristic fields. The resulting 
scheme was proven TVB for linear systems with both initial and initial-boundary con- 
ditions. Numerical results included both standard shock tube problems and a problem 
involving the interactions between a Mach 3 shock and a density wave - a prototype 
for shock-turbulence interactions. In all cases, we obtained results comparable to those 
obtained by recent nonoscillatory finite difference methods (e.g., [W], [S2]). The crucial 
generalization to multi-space dimensions was carried out in [C4], where we introduced a 
local projection limiting which does not have a direct counter-part in the current finite 
difference schemes, guarantees a local maximum principle for a class of very general trian- 
gulation (we called them B-triangulations), and maintains uniform high-order accuracy in 
the smooth part of the solution. Numerical results showed the potential of the scheme in 
easy handling general triangulations and boundary conditions. 

The organization of this paper is as follows. In §2 we briefly describe the formulation 
of the scheme. Special attention is paid to the monotone fluxes (approximate Riemann 
solvers) across the edges of the triangles and to the local projections associated with them, 
since they are distinct from the scalar case considered in [C4] . In §3 we present numerical 
results of our P 1 -RKDG scheme on the the shock reflection problem, the forward facing 


1 


step problem, and the double Mach reflection problem. We end with some concluding 
remarks in §4. 

2. The formulation of the scheme. Suppose we are solving the equation (1.1) 
( d = 2) on a polygonal domain £1. Let T/, = (K) be a B-triangulation of £1 (see [C4] for 
the definition), and set Vj, = {p 6 L°°(fl) : p\k G P k (K) VK e Tj} where P*(iiT) is the 
space of polynomials of total degree less than or equal to k on K. In this paper we only 
consider k = 1, which yields second-order accurate schemes. Notice that the functions in 
Vh can be discontinuous across the edges of the triangles in 7 *. For scalar equations, the 
discontinuous Galerkin method consists in finding u h : [0, T] i — * Vh satisfying 


( 2 . 1 ) 


/ ui l (t,x)vi i (x)dx + X / f( u h(t,x)) -n etK v h (x)dT 
K e£dK e 

- J f(uh(t,x)) • grad v h {x)dx = 0, Vvh G Vh, 

K 


where n ei x is the outward unit normal to the edge e and f = (/i,/ 2 )- This is obtained 
by multiplying (1.1) with a test function Vh G Vh, integrating over a triangle K € 7 *, 
integrating formally by parts, and replacing u by its approximation Uh- Since «/, can be 
discontinuous across an edge e, we replace f (uh(t,x)) • n e> if by a monotone flux 

(2.2) A.,*r(u»(l i "' (Rr) ), u h (x'‘« K >)), 


which is consistent: h e< K{u,u ) = f (u) • n e< K\ monotone: h e> K{u,v) is nondecreasing in u 
and nonincreasing in v; Lipschitz continuous; and legitimate as a flux: 

h'M'‘h(x in,(K) ),ui>(x”‘ ( - K) )) = -'!.,K'(f),(x i "‘ (K ' ) ),«l(i" ,(K ' ) )) for K'nK = c. 

The line integrals in (2.1) are replaced by the two-point Gauss quadrature rule, and the 
surface integrals replaced by the mid-point quadrature rule. The resulting O.D.E. is then 
discretized by the second-order accurate TVD Runge-Kutta method in [Si] and coupled 
with a local projection applied at the end of each inner step. This local projection is 
crucial for keeping the nonlinear stability (a local maximum principle) of the scheme; it is 
described in detail in [C4]. -- 

For systems of equations, (2.1) is satisfied by each component of u^. The monotone 
flux (2.2), however, is replaced by an approximate Riemann solver: the relevant quantities 
f (u h (x ,nt( - K \t) ■ n e<K and f(u h (x ext ( !< \t) ■ n Ci X axe left-multiplied by the left eigenvectors 
of some average Jacobian across e (we use the Roe average [Roe] associated with u k ar >d 
u K ,, the cell averages of u in K and K' whose common edge is e, in the normal direction 
to e), the scalar monotone flux (2.2) (we use the local Lax-Friedrichs flux [C2]) is then 
applied in each of the resulting local characteristic fields, and the results are projected back 
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by using the right eigenvectors of the same Jacobian. This is essentially a one-dimensional 
characteristics procedure across the edge e. Details about it can be found in [C3]. 

Next we describe the local projection. For the scalar case, a local maximum principle 
can be proven if we compute local approximate gradients using u/^ and u/c<, where K and 
K' are the immediate and some further away neighbors of K, and use them to limit the 
deviations of u h evaluated at Gaussian points from the mean uk-. The limiting is performed 
in such a way as not to destroy the formal accuracy in smooth regions. The details and 
proofs can be found in [C4] . For systems of equations, the limiting should be performed in 
the local characteristic fields. We adopt the simplified version of the local projection which 
involves only the midpoints of the edges of K instead of all the Gaussian points (as in [04] ) 
and limits the deviation by considering only the approximate gradient obtained by using 
u k, and u/c' 2 , where K'\ and K ' 2 are the two ‘forward’ facing triangles adjacent to 

iv; notice that the vector M - B k, where B# is the barycenter of the triangle K, is a 
positive combination of the vectors B k\ — B k and Bk"' 2 — B k, see Fig. 1. (In [C4] we 
used a projection that uses an additional approximate gradient obtained with backward 
facing elements K’\ and K' 2 . Our numerical experience shows that such a projection is 
not needed in this framework.) The relaxation factor b , which is the ratio allowed of the 
deviation of u h to ua' over that computed by the approximate gradient, is chosen to be 1.5 
( b > 1 guarantees second-order accuracy in smooth monotone regions). This simplified 
projection reduces the computational cost and seems to work well numerically. We again 
use the Roe average associated with u/^ and Uft-q, where KC\K'i = e, along the normal of 
e, to perform the characteristic decomposition for the projection limiting at the midpoint 
of e; see Fig. 1. The projection limiting is done at each midpoint independently form the 
projections at the other midpoints. Then, a simple readjustment is performed to enforce 
the conservativity of the method. This renders the projection limiting simple and efficient. 
In all our computations we took M = 50; see [C4]. 

3. The numerical examples. We use three standard test problems: the regular 
shock reflection, flow past a forward facing step, and the double Mach reflection problem, 
to display the behavior of our P J -RKDG scheme. We run our code on triangulations 
that are obtained, essentially, by taking the usual finite difference grid and dividing its 
rectangles into two triangles. (In principle, this should produce a bigger numerical diffusion 
and a stronger mesh distorsion; on the other hand, these undesirable effects could be 
counterbalanced by the inherent increase of degrees of freedom. Indeed, this seems to 
be the case.) The objective is to verify that our scheme gives results comparable to 
the results given by the recent nonoscillatory finite difference schemes for these standard 
test problems. Except for the second run in Example 1, we deliberately do not use the 
alignments of triangles with shock structures in order to carry out a fair comparison with 
finite difference schemes. For problems involving complicated geometries and/or boundary 
conditions, the advantages of our finite element schemes over finite difference methods 
would be more significant. Also, for problems involving both shocks and complicated flow 
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Fig. 1. Projecting u^k (M) at the midpoint M. The deviation Uh\x (M.)-uk 
is compared with b times the deviation at M obtained by the approximate gradient 
formed with U/c, U/c'i> and Ux* 2 . 

structures, P fc -RKDG schemes would show their advantage in high resolutions. Research 
along these directions is currently being carried out. 

In all our computations we have used the CFL-condition At sup |e|/[RT| < 0.3, 
where Xk is the modulus of the velocity plus the sound speed evaluated at the barycenter 
of the traingle K. We use the CRAY2 at the Minnesota Supercomputer Center to carry 
out our computations. The code is carefully written so that most of the major loops are 
fully vectorized. Due to the local structure of the algorithm, a parallel version of the 
CRAY2 code could be written; we plan to do this in the future. Our graphics has been 
done with the finite element package MODULEF. We have used 30 isovalues in all our 
contour figures. 

Example 1: Regular shock reflection. This well-known example involves an 
oblique shock reflecting from the lower boundary of a rectangular domain. The com- 
putational domain is 0 < x < 4.12829,0 <?/<!. The initial condition is p = \.q x = 
2.9, q y = 0,p = 1/1.4 throughout the computational domain. The boundary conditions 
applied are: inflow boundary condition at x = 0 (prescribe all components of u/, with 
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values same as the initial condition); outflow boundary condition at x = 4.12829 (all com- 
ponents axe freely flowing out); enforcing the post-shock condition at the upper boundary 
y = l(p = 1.7, q x = 2.61932, = -.506339, p = 1.52824); and enforcing a reflection 

boundary condition at the lower boundary y = 0 (if an edge e of a triangle is at the lower 
boundary y = 0, then the boundary value of uj, at (x, y = 0) is chosen to be the same 
as u h from y i— ► 0 + for p,q x , and p but opposite in sign for q y . The exact solution to 
this problem is an incoming shock of 29 degrees with the lower boundary and a reflected 
shock of 23.28 degrees. The exact solution past the second shock should be p = 2.68732, 
q x — 2.40148, q y = 0 and p = 2.93413. 

In Fig. 2b we show the pressure contour computed with a refinement of the triangu- 
lation in Fig. 2a corresponding to the uniform Cartesian grid ‘Ax = Ay = 1/40.’ Notice 
how the incident shock is approximated much better that the reflected one. This is due 
to the fact that the triangles are partially aligned with the incident shock and cut the 
reflected shock, increasing in this case the numerical diffusion; see Fig. 2b. In order to 
see how dramatically the result can improve when triangles are perfectly aligned with the 
shocks, we show in Fig. 2d the pressure contour computed with the triangulation in Fig. 
2c. In this case, the L 1 -error (with only three triangles) is only ten times bigger than the 
L 1 -error produced by the approximation displayed on Fig. 2.b (which uses 800 triangles). 

This test problem is simple in structure: three constant values separated by two shocks. 
We use it to test (i) the nonoscillatory property of our simplified local projection, (ii) the 
behavior of the numerical reflecting boundary condition at the lower boundary, and (iii) 
the effect when triangles are aligned with shocks. 



Fig. 2. a. The triangulation ‘Ax = Ay = 1/10’. The pressure below has been 
computed on a triangulation four times finer than this one. 
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Fig. 2.B. Pressure contour computed, on the triangulation ‘Ax = Ay = 



Fig. 2.C. The triangles of this iriangulation are aligned along the shocks of 
the exact solution. 



Fig. 2.D. Pressure contours computed on the triangulation of Fig. 2.c. 
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Example 2: Flow past a forward facing step. This is one of the two-dimensional 
problems Woodward and Colella [W] used to test the behavior of various finite difference 

f 0-6, if y < .2 

schemes. The computational domain is0<y<l,0<x<< . The 

l 3.0, ify > 0.2 

initial condition is a Mach 3 uniform flow: p = 1.4, q x — 3, q y — 0,p = 1 throughout the 
computational domain. The boundary conditions applied are: inflow at x = 0; outflow 
at x = 3; and reflecting at the walls y = 1, y = 0, x = 0.6 and y — 0.2. The corner 
is a singularity. In [W], Woodward and Colella suggested a way to numerically treat the 
singularity. We display results without, Fig. 3.b, and with, Fig. 3.c, such a treatment. 
Figs. 3.b and 3.c show the density contours at T — 4 computed with a refinement of the 
triangulation shown in Fig. 3. a corresponding to the Cartesian grid ‘Ax = Ay = 1/40’. 
We can see that our P^RKDG scheme produces results comparable with the same order 
MUSCL finite difference scheme using the grid ‘Ax = Ay = 1/40’.. We remark that since no 
sharpening technique is applied in the linearly degenerate field, the contact discontinuities 
are more seriously smeared than shocks. We are currently investigating the application of 
Yang’s artificial compression ideas [Y]. 

Example 3: Double Mach reflection. This is the second example used in [W] to 
compare various finite difference schemes. We use a different computational domain (Fig. 
4a) from the one used in [W]. Our domain is physically more natural and computationally 
easily manageable for triangulations. This is in contrast with finite difference methods 
which would meet complications in using our domain. The initial condition is described 
in Fig. 4a. It corresponds to a Mach 10 shock making 60 degrees with the bottom wall. 
The boundary conditions applied are: inflow at x = 0; outflow at x = 3; reflecting at the 
bottom; and the exact solution of a Mach 10 shock at the top. In Fig. 4c we show the 
density contour at T = 0.2. The tri angulation used is drawn in Fig. 4b. Again , we do not 
try to align the triangles with the shocks, and the number of degrees of freedom used is 
close to the middle case (Ax = Ay = 1/60) in [W]. We again observe a result comparable 
with the finite difference schemes, except for the above mentioned smearing of contact 
discontinuities. 

4. Concluding remarks. We have discussed the generalization of the RKDG method 
to two-dimensional systems of conservation laws, using the Euler equations of gas dynamics 
as an example. Numerical tests using P 1 elements on three standard examples show 
comparable results with nonoscillatory finite difference schemes and indicate good potential 
of our finite element method in handling geometry and boundary conditions. 
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Fig. 3. a. The triangulation ‘Ax = Ay = 1/5’. The densities below hve been 
computed on a triangulation eight times finer . 



Fig. 3.b. Density contours computed on the triangulation ‘Ax = Ay = 1/40 \ 
No treatment of the singularity at the corner is used . 



FlG. 3.C. Density contours computed on the triangulation ‘ Ax = Ay = 1/40’. 
The treatment of the singularity at the comer suggested in [W ] is used . 
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Fig. 4. a. The computational domain and the initial conditions. 




§&BB BBBB BBBBBBBBBBBBBBBBBBBBBB 
^ IBBBBB BBBBBBBBBBBBBBBBBBBBBBB 
^BBBBBBBBBBBBBBBBBBBBBBBBBBBB 




Fig. 4.b. The triangulation ‘Ai = Ay = 1/10’. The density below has been 
computed of a triangulation six times finer. 
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